function [rmax_coher] = Rmax_coherent_comb(iota, psi, phis, thetas, detNames, indigoDetPsi, mTot)


   f = linspace(0.1, 10000, 40000);
   f0 = min(f);
   logLambda = 8.0;
   m1 = mTot/2; % Equal mass binaries (for now)
   m2 = mTot/2; % 

   % Order of j : HLVI

   %===============================================================
   j = 1; % (Hanford)
   Sh = detnoisepsd('LIGO', f);
   inspDist = calcinspiraldist( m1, m2, f, Sh, f0); % in Mpc
   %The factor of 4 takes care of optimally oriented binaries (iota = 0.0)
   inspDist = inspDist/2;
   Det = LoadDetectorData (detNames(j));
   [Fp Fc] = ComputeAntennaResponse(phis, thetas, psi, Det.d);
   r1Sq = (inspDist/logLambda)^2 * ((1 + (cos(iota))^2)^2 * Fp.^2 + 4*(cos(iota))^2 * Fc.^2);
   
   %===============================================================
   j = 2; % (Livingston)
   Sh = detnoisepsd('LIGO', f);
   inspDist = calcinspiraldist( m1, m2, f, Sh, f0); % in Mpc
   %The factor of 4 takes care of optimally oriented binaries (iota = 0.0)
   inspDist = inspDist/2;
   Det = LoadDetectorData (detNames(j));
   [Fp Fc] = ComputeAntennaResponse(phis, thetas, psi, Det.d);
   r2Sq = (inspDist/logLambda)^2 * ((1 + (cos(iota))^2)^2 * Fp.^2 + 4*(cos(iota))^2 * Fc.^2);
            
   %===============================================================
   j = 3; % (Virgo)
   Sh = detnoisepsd('Virgo', f);
   inspDist = calcinspiraldist( m1, m2, f, Sh, f0); % in Mpc
   %The factor of 4 takes care of optimally oriented binaries (iota = 0.0)
   inspDist = inspDist/2;
   Det = LoadDetectorData (detNames(j));
   [Fp Fc] = ComputeAntennaResponse(phis, thetas, psi, Det.d);
   r3Sq = (inspDist/logLambda)^2 * ((1 + (cos(iota))^2)^2 * Fp.^2 + 4*(cos(iota))^2 * Fc.^2);

   %===============================================================
   j = 4; % (INDIGO)
   Sh = detnoisepsd('LIGO', f);
   inspDist = calcinspiraldist( m1, m2, f, Sh, f0); % in Mpc
   %The factor of 4 takes care of optimally oriented binaries (iota = 0.0)
   inspDist = inspDist/2;
   Det = LoadDetectorData (detNames(j));
   
   % This piece of code takes into account the arm geometry of the Indigo detector
      Det.psi = [indigoDetPsi 00 00]*[1; 1/60; 1/3600];

      theta = (90-Det.phi)*pi/180;
      phi   = (Det.lambda)*pi/180;
      psi   = (Det.psi*pi)/180;

      thetaCap = ...
              [ -cos(theta)*cos(phi); -cos(theta)*sin(phi); sin(theta) ];
      phiCap = [ -sin(phi); cos(phi); 0.0 ];

      Det.X    =  cos(psi)*thetaCap + sin(psi)*phiCap;
      Det.Y    = -sin(psi)*thetaCap + cos(psi)*phiCap;
      Det.Z    =  cross(Det.X,Det.Y);

      %----- Response matrix.
      Det.d = 0.5*(Det.X*Det.X' - Det.Y*Det.Y') ;

   [Fp Fc] = ComputeAntennaResponse(phis, thetas, psi, Det.d);
   r4Sq = (inspDist/logLambda)^2 * ((1 + (cos(iota))^2)^2 * Fp.^2 + 4*(cos(iota))^2 * Fc.^2);





  % Calculate the coherent rmax for 3IFO and 4IFO combinations 
   
   rmax_coher.HLV = sqrt( r1Sq + r2Sq + r3Sq );
   
   rmax_coher.HLVI = sqrt( r1Sq + r2Sq + r3Sq + r4Sq );
   
 
